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Spatial clustering of interacting bugs: 
Levy flights versus Gaussian jumps 
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Abstract. - A biological competition model where the individuals of the same species perform 
a two-dimensional Markovian continuous-time random walk and undergo reproduction and death 
is studied. The competition is introduced through the assumption that the reproduction rate de- 
pends on the crowding in the neighborhood. The spatial dynamics corresponds either to normal 
diffusion characterized by Gaussian jumps or to superdiffusion characterized by Levy flights, ft 
is observed that in both cases periodic patterns occur for appropriate parameters of the model, 
indicating that the general macroscopic collective behavior of the system is more strongly influ- 
enced by the competition for the resources than by the type of spatial dynamics. However, some 
differences arise that are discussed. 
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Introduction. — Interacting particle systems help to 
model and understand various problems in many diverse 
fields. In biological contexts they are particularly impor- 
tant to study aggregation phenomena of individuals. Fish 
schools, insect swarms, bacterial patterns, bird flocks and 
patchy plankton structures are just a few examples reveal- 
ing the ubiquity and fundamental importance of organism 
aggregates. 

An attempt to address a simple mechanism giving rise 
to the clustering of particles (with emphasis on plankton 
patchiness) was made within a Brownian bug model [T] 
(see also refs. [21I3])- The model consists of an ensemble 
'of particles (bugs), each one dying or reproducing with 
a given probability and undergoing Brownian motion. If 
the diffusivity of the particles is low enough, macroscopic 
spatial clustering occurs, since a newborn is located close 
to the parent but particles can die anywhere. If diffusivity 
is large, particles perform more extended walks and the 
region that was left empty due to the death of a particle, is 
occupied fast by some other particle. Similar results were 
obtained in refs. 0HZ] where lattice-models were studied. 

The basic Brownian bug model lacks any interaction 
between the particles. In a more realistic model (inter- 
acting Brownian bug model) the inter-particle interaction 
was taken into account assuming that the birth and death 
of individuals depend on the number of other bugs in 



the neighborhood [8H13]. For appropriate parameters, a 
salient property of that model is the formation of spa- 
tially periodic clustering of bugs [5]. For large diffusion 
the clusters become blurred and the periodic pattern is 
replaced by a more uniform distribution of bugs. Impor- 
tantly, whereas the positive correlations leading to clus- 
tering in the non-interacting bug case arise from the re- 
productive correlations, i.e., from the fact that offspring 
is born at the same location of parent [5,6, the periodic 
arrangement of clusters in the system of interacting bugs 
is a consequence of the competitive interaction and has a 
spatial scale determined by the interaction range R [8]. 

At the same time it has been observed that many liv- 
ing organisms move consistently with Levy flight behav- 
ior [14Ml8| . In particular, the motion of some bacteria is 
found to be described by Levy statistics [191I20]. as well 
as the movement of spider monkeys in search of food |21j . 
The Levy type of motion has been shown to be advanta- 
geous with respect to standard Brownian motion in some 
searching strategies involving foraging |18j . or in order to 
enhance encounter rates at low densities [35]. The main 
reason for this resides in the occurrence of occasional long 
jumps. 

However, the impact of Levy- type diffusion on the prop- 
erties of organism aggregates has not received much atten- 
tion thus far. In the present paper we investigate in the 
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context of a simple interacting bug model, continuous in 
space, how the occurrence of long jumps in the motion 
of the individuals influences the collective behavior (c.f. 
ref. [23 ). Similarities and differences between the inter- 
acting Brownian and Levy bug systems will be highlighted. 

Model and numerical algorithm. We consider a 
system consisting initially of Nq pointlike particles, placed 
randomly in a two-dimensional L x L square domain with 
periodic boundary conditions. After the random time r, 
a particle i, chosen randomly among all the N(t) bugs in 
the system at the present time t, undergoes one of the two 
events: it either reproduces or disappears. For the birth 
and death rates of the i-th particle we assume [8] , 



r\ = max (0, r bQ - aN l R ) 
r d = max (0, r d0 + PN R ) 



(1) 
(2) 



Here N R is the number of particles which are at a distance 
smaller than R from particle i, the parameters r^g and r d0 
are the zero-density birth and death rates, and the pa- 
rameters a and /3 determine how r\ and r d depend on the 
neighborhood; the function max() enforces the positivity 
of the rates. Such a choice for the reproduction and death 
rates introduces an interaction between the bugs. For the 
random times r an exponential probability density with 
the complementary cumulative distribution 



p(r) = exp(— r/f) 



(3) 



and a characteristic time (r) = f = i? to J is chosen; f is 
the time-scale parameter and 
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(rt+r d ) 



(4) 



is the total rate of all the N = N(t) particles. In the 
case of reproduction, the new bug is located at the same 
position (xi,yi) as the parent particle i. After the de- 
mographic event, i.e., after each random time r, all the 
particles perform a jump of random length £ in a random 
direction characterized by an angle uniformly distributed 
on (0, 2ir) (£ and the direction of the jump are different 
for each particle). The new present time is t + r and the 
process is repeated again, until the final simulation time 
is reached [13] . 

In order to simulate the system where the particles un- 
dergo normal diffusion, a Gaussian jump-length probabil- 
ity density function is used, 



<p{t) = (i 



<V2tt exp -t/{2l 2 ) 



(5) 



with variance (£ 2 ) — £ 2 ; £ is the space-scale parameter. 
Since we draw the angle specifying the direction of the 
jump from the interval (0,2-7r), we can neglect the sign 
of £. Note that the random walk defined in this way is 
not exactly the same as the one in which the walker per- 
forms jumps extracted from a two-dimensional Gaussian 



distribution, but it also leads to normal diffusion and al- 
lows a more direct comparison with the Levy case. The 
corresponding diffusion coefficient is 



K = (f)/(2(r}), 



(6) 



which expressed through the space- and time-scale pa- 
rameters reads, k = £ 2 /(2f). As we choose to fix the 
value of n, then the space-scale parameter is determined 
by £ — \J2kt = y / 2n/R to t- Note that because the number 
of particles is changing in time also the quantity i?tot is 
changing in time. However, as in this paper we investigate 
asymptotic statistically steady states, the number of par- 
ticles weakly fluctuates around a well-defined mean value. 
Thus, the space-scale parameter £ is fluctuating in time, 
but this time dependence is not affecting qualitatively the 
dynamics of the system. 

In order to simulate the system where the particles un- 
dergo superdiffusion one can use a symmetric Levy stable 
probability density function for the jump size, which be- 
haves asymptotically as [LHH7] 



-¥ ±00 



(7) 



with the Levy index < /i < 2. For all Levy stable 
probability density functions with \i < 2 the variance di- 
verges, (£ 2 ) — oo, leading to the occurrence of extremely 
long jumps, and typical trajectories are self-similar, on all 
scales showing clusters of shorter jumps intersparsed by 
long excursions. For < a < fi < 2 fractional moments 
converge, (|^| a ) < oo. Therefore, for the Levy index in 
the range 1 < /x < 2 the value of (|^|) is finite, whereas 
for < fi < 1 it diverges. The complementary cumulative 
distribution corresponding to (|7|) is 



£^±c 



(8) 



As a simple form of complementary cumulative distribu- 
tion function, which is normalizable and behaves asymp- 
totically as ©, we use 



P^£) = {l + b l i»\£\/£)-» , 



(9) 



where b = [T(l - fi/2)T(jj,/2)]/T(ji). As before, we can 
neglect the sign of £ since the direction of the jump is 
assigned by drawing an angle on (0, 2tt). Now the diffusion 
coefficient ([6]) is infinite, but one can define a generalized 
diffusion coefficient in terms of the scales £ and f as [14 , 17 



k„ = r/(2r) 



(10) 



Therefore, in the case of the Levy flights, when fixing the 
value of Kfj,, the space-scale parameter is, £ = (2k m t ) 1 / M = 
(2^/fitot) 1 ^. 

The model described can be interpreted in the follow- 
ing way: during a time interval r each individual moves 
on average a distance £ in some direction but only one 
of them reproduces or dies. For positive values of a and 
f3, the more neighbors a particle has within the radius 
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Fig. 1: Interacting Brownian bugs (a) versus interacting Levy 
bugs (b): spatial configuration of particles in the statistically 
steady state. The parameters are: rj,o = 1, rdo = 0.1, a = 0.02, 
ft = 0, and R = 0.1 (see also text). In the case of Brownian 
bugs the diffusion coefficient is k = 10 -5 . In the case of Levy 
bugs the value of Levy index is fj. = 1 and the anomalous 
diffusion coefficient is k u = 56 x 10" 5 . The average number of 
particles in the two systems is approximately the same, (TV) = 
2555 and (TV) = 2565, respectively. 



R, the smaller is the probability of reproduction and the 
larger is the probability that the bug does not survive, 
e.g., due to competition for resources. The periodicity of 
the simulation domain represents the fact that particles 
can leave and enter the observation area. In the case of 
long jumps it represents the situation where there is one 
particle that leaves the domain and another one that ar- 
rives from somewhere, perhaps from far away, and its po- 
sition is basically random in the domain. In principle, one 
could also simulate the system so that after the inter-event 
time r a randomly chosen particle undergoes one of the 
three events: reproduction, death, or jump. However, the 
results obtained are the same as following the procedure 
described above; one would just have one more parameter, 
the jump rate, and different numerical values for the other 
parameters. 



Results. — In the following we set (3 — in eq. © 
and a = 0.02 in eq. ((T|), i.e., the probability of death is 
constant and the same for all the particles while the prob- 
ability for the reproduction depends on how crowded the 
environment is. Also, we use the values r&o = 1, ^do =0.1, 
L = 1, and R = 0.1 in all the figures presented in the 
current paper. For ft = the critical number of neigh- 
bors, N R , for which death and reproduction are equally 
probable for particle i, is determined by 



N R = A M /a , 



(11) 



where Am = r^o — r c i$. If N R < N R it is more probable 
that particle i reproduces and if N l R > N R death is more 
likely. In the following we fix the value of the Levy index 
to jU = 1. The detailed discussion of the influence of the 
birth and death rates as well as the Levy index on the 
results is outside the scope of the current paper and will 
be presented elsewhere. 

In the case of interacting Brownian bugs, for small 
enough k and large enough A^, the occurrence of pe- 
riodic patterns has been observed previously [3[T0]. In 
the statistically steady state clusters form and arrange 
in a hexagonal lattice (see fig. [TJa). For large values of 
the diffusion coefficient the periodic pattern is replaced 
by an almost homogeneous distribution of particles (see 
also fig. [5]-a). In the case of Levy flights, since the dif- 
fusion coefficient ([6]) is infinite, one could expect that for 
the interacting Levy bugs the spatial distribution will not 
reveal a periodic pattern. However, as can be seen from 
fig. [TJb, this is not the case. The reason for the diver- 
gence of the diffusion coefficient in the Levy case is in the 
statistical weight of the large jumps. These large jumps 
have some influence in the characteristics of the patterns 
formed, but the large-scale structure is ruled mainly by 
the interactions between particles. 

For the interacting bug system with Gaussian jumps and 
moderate diffusion the appearance of periodic clustering 
is well captured by combining the effects of diffusion and 
interaction in a mean- field approach j8]. Following the 
steps in ref. [8] one obtains the following equation as a 
mean-field approximation to the dynamics of the density of 
particles p(x, t) in the interacting Levy bug model (j3 = 0): 



flp(x,t) 
dt 



p(x,t) (am-v J dyp(y,t)^j +K M V 7 p(x,t) . 

(12) 

The integration domain D is the set of points within a dis- 
tance smaller than R from x: |x — y| < R. The operator 
V 7 , with 7 = 43, is the fractional diffusion operator asso- 
ciated to the Levy flights of exponent fi G (0, 2) [141IT7] . 
It reduces to the standard diffusion operator for fj, > 2. In 
this mean-field description the first term accounts for the 
net growth of the population, the second one takes into 



1 This expression was incorrect in the paper originally published 
in EPL 92, 40011 (2010); it was corrected to the present form in the 
Erratum in EPL 95, 69902 (2011). 
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Fig. 2: Comparison of the radial distribution functions for 
the systems with Gaussian jumps and with Levy flights. The 
parameters are as in fig. [T] 



account the non-local contribution associated to the sat- 
uration due to the interactions within a distance R, and 
the third term describes the spatial diffusion of particles. 
Equation (|12l) differs from the mean-field approximation 
derived in ref. [5] only in the third term where the stan- 
dard diffusion operator is replaced by the fractional diffu- 
sion operator and the diffusion coefficient by the general- 
ized diffusion coefficient. Similar descriptions of reaction- 
diffusion systems with a fractional diffusion term can be 
consulted, e.g., in refs. [I4 l lll p 4 p 5] . 

Equation ([T2|) has always the uniform solution 



p(x,t) = pa = A M /(cnri? 2 ) , 



(13) 



which turns out to become unstable for small and/or 
large (what is small or large depends also on the value 
of p). This can be seen by introducing the ansatz p(x, t) = 
po + Sp(x,t) in eq. (|T2|) and linearizing in Sp(x,t). The 
result is 

<S/j(x, t) ~ exp(ik • x + A(k)t) , (14) 



with 



A(|k|) = -« M |k|T - 2A M J 1 (|k|i?)/(|k|i?) ; 



(15) 



J\ is a Bessel function. An instability of the uniform solu- 
tion will occur if the sign of A(|k|) changes from negative 
to positive at some value of |k|. This will occur first for 
the critical wavenumber |k c | for which A(|k|) is maximum. 
The instability will develop in a periodic pattern which, at 
least close enough to the instability, will have a periodic- 
ity (5 = 27r/|k c |. Introducing a dimensionless wavenumber 
q = |k|i? and growth rate A = WX/k^, eq. (ITBl reads, 



A(9) = -g* -vJi(q)/q 



(16) 



The latter equation shows that, within this mean-field de- 
scription, the relevant parameters are 7 (or p) and the 
dimensionless quantity v = 2R 1 A a {,/ n^. Imposing the 
condition A(q) — 0, corresponding to the change of sign in 
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Fig. 3: Radial distribution function for the system with (a) 
Gaussian jumps and (b) Levy jumps for several values of the 
diffusion and generalized diffusion coefficient, respectively. The 
other parameters are as in fig. [J] (see also text). 



the growth rate, and A'(q) — 0, corresponding to the max- 
imum of the growth rate, one can find the critical value 
of v, v c — —q2 +1 1 ' J\{qc), and the equation for the critical 
wavenumber: 



qcJ2{q c )/Ji{q c ) = -7 ■ 



(17) 



Numerical solution of the latter equation provides q c , and 
therefore 5, as a function of 7 (or p), and that 5 is propor- 
tional to R. In order to make a comparison with the results 
from numerical simulations of the system, one should keep 
in mind that the periodic boundary conditions will change 
q c to the closest number of the form (2nR/L)\/n 2 + m 2 
with n and m integers. For p = 1 the value of 8 ob- 
tained from the mean-held approximation is 0.128036 and 
for p = 2 (Gaussian case) it is 0.131306, i.e., the periodic- 
ity of the pattern is in both cases of the order of R = 0.1. 

To extract the periodicity of the pattern from the nu- 
merical simulations of the system, it is useful to study 
the radial distribution function g(r). It describes how 
the density varies with the distance from a given parti- 
cle respect to the one expected from a uniform distribu- 
tion, giving thus additional information about the distri- 
bution of bugs [26] . It is computed in the standard way, 
i.e., by counting all particles, dn, at a distance between 
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Fig. 4: The distribution of cluster sizes for the interacting 
Brownian bug and Levy bug models. All parameters have the 
same values as in figs. [1] and [2] 

r and r + dr from the target particle, using the formula 
dn = (N / L 2 )g(r)2irrdr , and averaging over all particles 
and different long times. In fig.[2]the comparison of the ra- 
dial distribution functions for the Brownian and Levy bugs 
(for the same systems as in fig. [TJ is depicted. Figure [3}a 
shows the behavior of g(r) in the case of Brownian bugs 
and in fig. [3]L in the case of Levy bugs for various values 
of diffusion coefficient and generalized diffusion coefficient, 
respectively. The second maximum of g(r), indicating the 
periodicity of the pattern, appears in fig. [3}a (in the sys- 
tems of Brownian bugs), up to k = 10 -5 , at 0.13125 and in 
fig-Gl-b (in the Levy case), up to k p = 10~ 3 , at 0.12875, be- 
ing in good agreement with the results obtained from the 
mean-fieled approximation, and becoming slightly larger 
only at larger values of k and k^. The anomalous expo- 
nent fi has only a very light influence on the periodicity 
of the pattern, which is of the order of R (see also fig. [2] 
where the slight shift can be can be noticed). 

Figure [3] shows also that as K or increases, the first 
peak of g(r) gets lower and wider: the total number of 
particles in the system decreases and the clusters become 
more spread. For large values of diffusivity (k or *t^), 
the oscillations of g(r) smooth out, indicating that the 
periodic pattern becomes replaced by a more homogeneous 
distribution of bugs. 

As a difference compared to the case of interacting 
Brownian bugs, we observe that now, even at small val- 
ues of Kfj,, there are many solitary particles appearing for 
short time periods in the space between the periodically 
arranged clusters, due to the large jumps, c.f. figs. [TJa 
and [JJb. This is even better illustrated by fig. 2] where 
the probability distributions of the cluster sizes for the 
systems with Gaussian and Levy jumps are depicted (for 
the same systems as in fig. [TJ i.e., for the given param- 
eters the values of k and have been chosen such that 
the average number of the particles in the two systems is 
approximately equal). The clusters are defined using the 
nearest neighbor clustering method [27] , with a threshold 
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Fig. 5: Cross-section of the two-dimensional particle density 
of the average cluster in (a) linear and (b) semi-log scale. Com- 
parison between the interacting Brownian and Levy bug mod- 
els. The values of the parameters are the same as in figs. [TJ 
[2l and [4] The continuous lines correspond to the fitting with 
Gaussian functions. 



value 0.025 to define neighbors. As one can see, for the 
system of interacting Brownian bugs the most probable 
size of large clusters forming the periodic pattern is 41 
particles whereas for interacting Levy bugs it is 36 (the 
critical number of neighbors is in both cases = 45). 
A noticeable difference occurs in the proportion of single- 
particle clusters: the peak of the distribution at the value 
characterizing isolated particles (N c = 1) is much higher 
in the Levy flight case. Returning to fig. [5J which refers 
to the same systems as in figs. [TJ and 2J it is to be noticed 
that the minimum between the first and the second peak is 
much lower for the system of Brownian bugs. This result 
is consistent with fig. [4] in the system with Levy flights, 
the occurrence of many single particles between the clus- 
ters forming the periodic pattern produces a higher inter- 
cluster density is observed. 

Finally, in fig.[5]the cross-section of the two-dimensional 
particle density of the average cluster is depicted for the in- 
teracting Brownian and Levy bug systems (same systems 
as in fig. [JJ) . The average cluster is obtained by setting 
the origin at the center of mass of each cluster forming the 
periodic pattern and averaging over all the clusters in the 
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simulation area and over time. In the case of Brownian as 
well as of Levy bugs, the central part of the average clus- 
ter, where most of the particles are concentrated, is well 
fitted by a Gaussian function (continuous curve). How- 
ever, fig. [5}b reveals the difference in the way the particle 
density decreases in the two systems when moving away 
from the center of mass of the cluster. Though even in 
the case of the Brownian bugs the tail of the average clus- 
ter is not really well described by a Gaussian function, a 
Gaussian decay provides a good approximation. Instead, 
in the case of the Levy bugs the tail of the average clus- 
ter decays slower. Since the average cluster is calculated 
from individual clusters which are separated by a distance 
of order R, one cannot properly estimate an asymptotic 
decay, but the decay in the Levy case seems close to ex- 
ponential. In any case, it is distinctively faster than the 
power-law behavior which would be exhibited by clusters 
of non-interacting particles moving purely by Levy flights. 
The existence of the interaction range R introduces a cut- 
off distance which makes the long-jumps characteristic of 
Levy flights not so relevant to determine the large-scale 
properties of the spatial particle configurations. 

Conclusion. — In the present paper we have studied 
a simple birth-death model with competition among the 
individuals of the same species. In particular, we have 
investigated how the superdiffusive motion of the individ- 
uals characterized by Levy flights influences the collective 
behavior. We have observed that the appearance of peri- 
odic clustering known from previous studies of the inter- 
acting Brownian bug model takes place also in the inter- 
acting Levy bug model for appropriate parameters. This 
is against the expectation that particles performing Levy 
flights cannot give rise to space-periodic clustering since 
their diffusion coefficient is infinite, for which the interact- 
ing Brownian bug model is known to reveal no periodic 
pattern. However, as a difference we have observed that, 
in the interacting Levy bug model, due to the long jumps, 
there are many single particles between the clusters, lead- 
ing also to the differences in the particle density profiles 
of the average cluster or in the short-distance behavior of 
the radial distribution function. From this one can con- 
clude that the large-scale collective behavior of the system 
is much more strongly influenced by the competitive in- 
teraction than by the type of spatial motion performed by 
the bugs. We have also verified that the mean- field ap- 
proximation (|I 2|) is proper to describe the periodicity in 
the interacting Levy bug model. 

As a final remark, let us mention that though the model 
studied in the present article describes rather living organ- 
isms, such as animals or bacteria, non-local interactions 
and Levy flights are important also in plant ecology for 
the development of vegetation patterns [28], due to the 
competition for the resources and because the seed dis- 
persal is often described better by Levy flights than by 
Gaussian jumps [29] . 
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